Exploring the effectiveness of flavone derivatives for treating liver diseases: Utilizing DFT, molecular docking, and molecular dynamics techniques

In exploring nature's potential in addressing liver-related conditions, this study investigates the therapeutic capabilities of flavonoids. Utilizing in silico methodologies, we focus on flavone and its analogs (1–14) to assess their therapeutic potential in treating liver diseases. Molecular change calculations using density functional theory (DFT) were conducted on these compounds, accompanied by an evaluation of each analog's physiochemical and biochemical properties. The study further assesses these flavonoids' binding effectiveness and locations through molecular docking studies against six target proteins associated with human cancer. Tropoflavin and taxifolin served as reference drugs. The structurally modified flavone analogs (1–14) displayed a broad range of binding affinities, ranging from -7.0 to -9.4 kcal mol⁻¹, surpassing the reference drugs. Notably, flavonoid (7) exhibited significantly higher binding affinities with proteins Nrf2 (PDB:1 × 2 J) and DCK (PDB:1 × 2 J) (-9.4 and -8.1 kcal mol⁻¹) compared to tropoflavin (-9.3 and -8.0 kcal mol⁻¹) and taxifolin (-9.4 and -7.1 kcal mol⁻¹), respectively. Molecular dynamics (MD) simulations revealed that the docked complexes had a root mean square deviation (RMSD) value ranging from 0.05 to 0.2 nm and a root mean square fluctuation (RMSF) value between 0.35 and 1.3 nm during perturbation. The study concludes that 5,7-dihydroxyflavone (7) shows substantial promise as a potential therapeutic agent for liver-related conditions. However, further validation through in vitro and in vivo studies is necessary. Key insights from this study include:• Screening of flavanols and their derivatives to determine pharmacological and bioactive properties using ADMET, molinspiration, and pass prediction analysis.• Docking of shortlisted flavone derivatives with proteins having essential functions.• Analysis of the best protein-flavonoid docked complexes using molecular dynamics simulation to determine the flavonoid's efficiency and stability within a system.

• Screening of flavanols and their derivatives to determine pharmacological and bioactive properties using ADMET, molinspiration, and pass prediction analysis.• Docking of shortlisted flavone derivatives with proteins having essential functions.
• Analysis of the best protein-flavonoid docked complexes using molecular dynamics simulation to determine the flavonoid's efficiency and stability within a system.novel synthetic derivatives of chalcone and flavone hybrids with 1, 2, and 3 triazole linkages for antimalarial activity found that optimization improved flavone binding affinity, resulting in higher antimalarial effects and more potent inhibition [35][36][37] .
In our research we investigated six major proteins, each selected based on specific criteria and their potential relevance to liverrelated conditions.The EGFR protein (PDB: 4UV7) was identified as the first protein of interest, inspired by Berasain and Avila's study, which suggested a potential role of EGFR in liver disease, particularly in conditions involving liver cell death.The link between EGFR and liver disease is supported by findings indicating that EGFR signaling provides protection against hepatocellular apoptosis induced by Fas ligand during acute liver injury and complications [38] .Furthermore, the HER2 protein (PDB: 7XJH) emerged as another focal point due to its potential involvement in liver disease, specifically hepatocellular carcinoma (HCC), the most prevalent form of primary liver cancer.HER2 overexpression has been observed in HCC, suggesting its potential as a therapeutic target for this disease.Shi JH and her team's work on HCC highlighted the correlation between HER2 overexpression and tumor stage, with 82% of samples exhibiting this relationship [39] .
In addition to EGFR and HER2, the FPPS protein (PDB: 1YQ7), a vital component of the mevalonate pathway crucial for cell growth, was included in our investigation.Alterations in FPPS could affect liver function or contribute to liver disease.Another enzyme, HPGDS (PDB: 1V40), responsible for producing PGD2, a mediator of inflammation, was considered due to its expression in the liver and potential involvement in inflammatory liver diseases, as indicated in the Human Protein Atlas.DCK (PDB: 1p60), an enzyme pivotal for activating nucleoside analogs in cancer and antiviral therapy, was also among our selected proteins.Its expression and activity may be altered in liver diseases, with cirrhosis, for instance, associated with lower DCK expression and sensitivity [40] .Lastly, the Keap1-Nrf2 signaling protein (PDB: 1 × 2 J) was chosen for its linkage to liver disease.This signaling pathway impacts various liver diseases, including hepatitis, cirrhosis, and liver cancer.Keap1-Nrf2 signaling controls genes related to detoxification, inflammation, and cell survival.Disruptions in this signaling pathway can lead to increased or decreased Nrf2 activity, with diverse impacts on liver disease development and outcome [41] .
Our study aimed to investigate flavone pharmacology, toxicity profiles, and biological activities and its 13 derivatives ( 1 ̶ 14 ).We employed density functional theory (DFT) to assess their stability and molecular properties, including thermodynamic stability, HOMO-LUMO gap, hardness, softness, chemical potential, electrophilicity index, and dipole moment.We also analyzed molecular docking to evaluate their binding affinities against six targeted proteins using 7,8-dihydroxyflavone and dihydroquercetin as the reference drugs.We assessed each flavone derivative to predict their activity spectra for substances (PASS).We investigated how these analogs interact with various key proteins in the body using molecular docking.We performed molecular simulations to validate our findings and assess the entropic ability of the drug candidates.Furthermore, we carried out an MD simulation analysis on the active site of the protein -ligand interaction to evaluate the stability of the protein -ligand complex.Our study aimed to develop a novel class of targeted anti-liver disease drugs with significant biological action.
The flowchart of the study is shown below:

Computational analysis
To explore potential variations in flavone and its derivatives, we initiated a screening process for compounds labeled ( 1 ̶ 14 ) within the PubChem molecular database.We then acquired their SDF structures and employed ChemDraw to depict these derivatives, as illustrated in Fig. 1 visually.This method facilitated the identification of isomers and their divergence, ultimately creating innovative and novel compounds.Subsequently, optimization procedures were conducted using Gaussian16 to refine the geometries of the structures ( 1 ̶ 14 ) [42] .These optimizations are detailed in Tables S1 to S14 and were performed at different levels of theory, specifically B3LYP/6-31G(d,p), B3LYP/6-311G(d,p), and B3LYP/6-311 ++ G (d,p) [43][44][45] .Frequency calculations were also carried out to confirm that all structures were in their energy minima.Previous studies [46][47][48][49][50][51][52] have demonstrated the accuracy and efficiency of using the B3LYP/6-31G(d,p) level of theory.The analysis of frontier molecular orbital (FMO) energy distribution, spanning from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO), along with the generation of molecular electrostatic potential (MEP) maps, was performed using GaussView 6 software.For a thorough visual presentation, the Supplementary Information (SI) encompasses optimized HOMO, LUMO, MEP maps, Natural Bond Orbital (NBO) data, and all structures pertinent to docking and molecular dynamics (MD), as displayed in Figs.S1 to S18 within the SI.
where μ and η were used to derive the electrophilicity indexes.Furthermore, to obtain accurate energy values and atomic charges, the single point energy and natural bond orbital (NBO) are used [57 , 58] .

Analysis of physicochemical and pharmacokinetic properties
To ensure that a drug achieves an appropriate concentration at its target site, ADMET properties and a drug's adherence to the Lipinski rule must be evaluated.In this study, three web-based tools, SwissADME, AdmetSAR, and PKCSM, were used to predict the ADMET properties of flavonoids and check if they violated any rule.SwissADME utilizes powerful predictive abilities to determine the physicochemical properties of drug candidates, while AdmetSAR employs five quantitative regression models to provide the most accurate drug prediction outcome.PKCSM uses a distance-based approach to compute a drug's potency, safety, and pharmacokinetic properties.The physicochemical and drug-likeness properties of the flavonoids are presented in Table 2 , and the canonical SMILES of each flavonoid were entered into both server tools to obtain the ADMET values.Additionally, the Molinspiration Chemoinformatics server tool was utilized to predict the downstream biological activity of flavonoids when interacting with various proteins in different target pathways, such as GPCR, KI, ICM, NRL, PI, and EI [58][59][60] .

Pharmacological properties
The Prediction of Activity Spectra of Substances (PASS) server tool can predict various pharmacological activities of compounds, such as toxicity, antibacterial, antifungal, mutagenicity, and more [61] This tool can predict approximately 300 pharmacological effects and can aid in identifying potential organic drug candidates.Pass Prediction was employed to analyze the flavonoids in this study, and the Pass Online Server tool was used for this purpose [62][63][64] This web tool provides access to biochemical information on compounds approved for medicinal use in the USA and Russian Federation [63 , 64] .The results of the analysis are presented in Table S15 in the SI.

Preparation of target protein
In order to assess the efficacy of a drug in preventing the progression of liver disease, it is crucial to determine its capacity to bind to a cancer-promoting receptor or oncoprotein.This study focused on six critical proteins, including epidermal growth factor protein (PDB: 4UV7), HER2 protein (PDB: 7JXH), farnesyl diphosphate synthase protein (PDB: 1YQ7), hematopoietic prostaglandin D synthase (PDB: 1V40), deoxycytidine kinase protein (PDB: 1P60), and Keap1 on Nrf2 repressor protein (PDB: 1 × 2 J).The 3D crystal structures of these proteins were obtained from the RCSB Protein Data Bank (PDB), the largest archive of thousands of sequenced proteins available for research in various formats [65] .The crystal structures were validated using X-ray diffraction to ensure the quality of the primary protein structure.The structures of all flavonoids were also optimized using UCSF Chimera software version 1.16 [66] .The histidine protonation state was set, and the standard residue was kept at the default AMBER ff14SB for both the polysaccharides and proteins.The Gasteiger method was used as the charged method and applied to other residues of each protein.
The protein structure was prepared by removing ligands, water molecules, and metal ions.Finally, Kollman charges were added to the protein molecule before being converted to PDBQT format by the AutoDock Tools (v.1.5.6) package of Pyrx software.The active sites of all the proteins were determined using the CastP server tool.The active sites of the two proteins are shown in Table 5 .The rest are given in the supplementary information.The PDB ID: 1p60 was chosen as it is a remarkable liver disease biomarker [67] .Keap1 on the Nrf2 repressor protein (PDB: 1 × 2 J) helps to elevate the accumulation of Nrf2 in the nucleus and protects hepatocytes against acute drug toxicity and inflammatory liver injuries [68] .Several liver diseases are associated with increased hepatic farnesyl diphosphate synthase in mice.Moreover, (PDB: IYQ7) is also suppressed when liver polyunsaturated fatty acids increase [69 , 70] (PDB: 1V40) is associated with hepatic inflammation and has therapeutic potential [71] (PDB: 4UV7) has shown increased activity and associations with liver disease therapy.Finally, recent studies have found that the Her2 protein (PDB: 7XHJ) has been shown to have a positive correlation with liver diseases [72] .

Protein-ligand docking
The docking simulations were carried out using PyRx version 0.8 software with AutoDock 4 and AutoDock Vina tools for ligand -protein interaction [68][69][70][71] .Prior to docking, the flavonoid structures were energy minimized with the Gasteiger method using UCSF Chimera software and then converted to PDBQT format using PyRx.For complete protein surface coverage during docking, a grid box was set up for each protein.The protein -ligand complexes were stable throughout the docking process and were verified by redocking.UCSF Chimera software was used to identify amino acid interactions at the docking site, while PyMOL version 2.5 was used to generate molecular docking images [72] .Moreover, BIOVIA Discovery studio was used to recognize 2D protein interactions between the ligand and protein and to examine the hydrogen density surrounding the interacting residues with the protein [73] .The results of the docking studies are presented in Table 4 .

Molecular dynamics simulation
To investigate the interactions between proteins and ligands and better understand the structure and function of biological macromolecules, molecular dynamics (MD) simulations were employed using the GROMACS version 2021.6 package [74][75][76] .The AM-BER99SB force field was used to describe atomic interactions in the system.Although computationally intensive, MD simulations provide more precise insight into complex systems than docking and can reveal features inaccessible by experimental methods [77 , 78] .The complex with the lowest binding and negative binding affinity for flavonoid ( 7) was used to conduct MD simulations for PDB: 1 × 2 J and PDB: 1P60.Protein topology parameters using the Galaxy European Server [79] .The system was solvated in a triclinic box of SPC water molecules and neutralized with Na + ions to achieve default salt concentrations [79][80][81] .Using the leap-frog algorithm, we equilibrated the system through a position-restrained dynamic simulation (NVT) at 300 K for 3000 ps.A production run was conducted for an additional 3000 ps at the same temperature and pressure settings.The MD trajectories were visualized and analyzed using VMD, PyMol and GROMACS programs [74][75][76][77][78][79][80][81] .MD simulation for the docked compounds was performed for 100 ns at 298 K.Moreover, gmx rms, gmx gyrate, and gmx rmsf utilities were used for root-mean-square-deviation (RMSD), radius-of-gyration (Rg), and root-mean-square fluctuation (RMSF) analysis, respectively.We also employed the GROMACS utility 'gmx hbond' to perform hydrogen bond analysis, and the temperature and potential energy were determined.Additionally, the Bio3D package through the Galaxy European Server was employed for the MD trajectories' principal component analysis (PCA) [77][78][79][80][81][82][83] .

Analysis of frontier molecular orbitals using DFT
The reactivity and stability properties of flavone and its analogs ( 1 ̶ 14 ) were assessed using FMO analysis, which provided valuable insights.FMO analysis helped determine the chemical reactivity, softness, hardness, chemical potential, and electrophilic index of the compounds through the HOMO-LUMO energy gap (Egap) parameter.A large Egap indicates high stability and low chemical reactivity, whereas a narrow Egap corresponds to high chemical reactivity and low stability.The energy levels of HOMO and LUMO also play a crucial role in understanding the compounds' electron acceptor and donor properties, respectively, and are critical in chemical reactions.Therefore, FMO analysis is essential for evaluating the reactivity and stability properties of molecules.Furthermore, the ionization potential (IP), electron affinity (EA), electronegativity ( ), chemical potential (μ), global hardness ( ), softness ( ), electrophilicity (  ), and dipole moment of all compounds ( 1 ̶ 14 ) were calculated using the B3LYP/6-31G(d,p) level and are presented in Table 1a (see Figs. S1 to S18 in the SI).The Hartree-Fock model describes the movement of electrons from the ground state to the excited state in atomic orbitals during chemical reactions in a biochemical system [84] .A more significant energy gap indicates a lower electronic transition of the electron, resulting in decreased chemical reactivity of the molecules.Conversely, a smaller gap exhibits a higher atomic system and chemical reactivity, closely related to binding with the target protein [84] .In Table 1, the energy gaps of all compounds fall within the range of 3.57 eV to 4.11 eV.Adding halogens to flavonoids did not result in any significant difference in the energy gap values compared to mono-and di-hydroxyflavones, indicating that group 7 elements do not play an essential role in the binding site of the receptor protein.The investigation of the compounds reveals that they all display a higher hardness ( ) than softness ( ).The electrophilicity index, an essential parameter for pinpointing reactive sites, depends on the electron transfer between the acceptor (LUMO) and the donor (HOMO), along with the corresponding energy changes [84][85][86] .Take flavonoid ( 7 ) as an example, its electrophilicity index is 3.946 eV.This measurement acts as a benchmark for evaluating the biological activity of flavonoids ( 7) and gauging their performance in molecular docking studies.Notably, flavonoid ( 7 ) has a higher eV value than the reference drug, which is 3.913 eV.It also shows a lower Debye value of 3.688, signifying its effectiveness.These findings indicate a stronger interaction between flavonoid ( 7 ) and the active site amino acids of the target protein, reinforcing its potential as a potent agent.
Structural parameters, encompassing bond lengths and angles for these compounds, were computed using the B3LYP level of theory with varying basis sets (6-31G(d,p), 6-311G(d,p), and 6-311 ++ G (d,p)).To assess the precision of the selected method, we conducted a comparative analysis between the computational results and the experimental data.This comparison aimed to identify discrepancies or variations between the calculated and experimental results [87][88][89] .Table 1b displays discrepancies in structural parameters resulting from methodological disparities in the analysis of single crystals of commercially available compounds, namely, compounds ( 1) and ( 7) .These crystals were obtained through room temperature X-ray diffraction, and detailed information on crystal data, data collection, and structure refinement can be found in previously published reports [87][88][89] .Notable differences  were observed when comparing the optimized structures with the crystal structures of these compounds.Specifically, variations were identified in the position of the carbon six (C5-C6 double bond) in the inner ring, as well as differing conformations in the outer ring, particularly in the C1-C2 and C1-C6 regions, along with bond angles, as indicated in Table 1b .These variations resulted in calculated bond lengths within the ring structure that deviated from the crystallographic data.The electrostatic potential maps offer qualitative insights into the net electrostatic impact at a point because of the total charge distribution.These maps help identify electrophilic and nucleophilic sites, where electrophilic sites act as electron acceptors and nucleophilic sites act as electron donors.ESP contour maps demonstrate the electron distribution that favors a low potential for high wavelengths and a high potential for low wavelengths.For all the compounds analyzed, regions of respective colors indicating electrophilic (blue) and nucleophilic (red) sites and partial nucleophilic sites (yellow regions) were identified, indicating total density [85] .The MEP maps of the molecules, including Fig. 2 and Figs.S12-S14 in the SI, can help determine the regions of electrophilicity and nucleophilicity.

Analysis of the pharmacological characteristics of flavonoids
In this study, multilevel neighborhoods of atoms (MNA) descriptors in the prediction of activity spectra for substances (PASS) were employed to conduct a comprehensive analysis of the pharmacological properties of 14 flavone derivatives ( 1 ̶ 14 ).Flavonoids are known for their diverse biological activities and ability to modulate reactions by binding to downstream proteins.The PASS predictions for each compound were based on various parameters, including their potential to affect membrane integrity, act as membrane permeability inhibitors, exhibit anticarcinogenic properties, and display antioxidant activities.The Pa and Pi values for these parameters ranged from 0.001 to 1.000.A Pa value exceeding 0.6 indicated significant biological activity in in vitro experiments.
Notably, flavonoids ( 1, 7 , and 9) demonstrated Pa values ranging from 0.618 to 0.989 for all four parameters, as detailed in Table S15 in the Supplementary Information (SI).These findings suggest the potential for these compounds to maintain good cell membrane integrity and block cell permeability effectively.Flavonoid ( 9) exhibited higher values than the control drug and exhibited noncarcinogenic and antioxidant properties.However, the membrane integrity and permeability values for Flavonoid ( 11) could not be determined.Furthermore, flavonoids ( 1, 2, 6, 8, 9 , and 11) scored below the cutoff value for anticarcinogenic activity, indicating that they may be susceptible to carcinogenic substances.Nonetheless, all compounds demonstrated scores above the cutoff value for being membrane integrity agonists and membrane permeability inhibitors, suggesting their effectiveness in these two parameters.In summary, the PASS Prediction analysis provides valuable insights into the potential biological activities of the studied flavonoids.This information can assist in identifying promising candidates for further investigation and drug development.Cell survival, structural integrity, and polarity heavily rely on the stability and functionality of cell membranes.Compounds that serve as membrane integrity agonists or contribute to membrane stabilization can avert apoptosis, metabolic complications, and nonalcoholic fatty acid disease.Recent research has identified the targeting of mitochondrial membrane permeability as a promising pharmacological approach for liver diseases, including fatty acid liver disease and cardiovascular disorders [90][91][92] .Besides these attributes, an effective drug should have anticarcinogenic and antioxidant properties to neutralize radical formation and aid chemotherapy.Flavonoids, especially those containing ortho-dihydroxy groups on the C ring of flavones, have shown significant effectiveness in scavenging radicals, making them crucial for antioxidant and anticancer properties [92 , 93] .In conclusion, compounds that stabilize cell membranes or function as membrane integrity agonists are essential for cell survival, and those with antioxidant and anticancer properties can be beneficial supplements to chemotherapy [92] .

Analysis of physiochemical and drug likeness properties of flavonoids
The physicochemical properties must be examined to assess whether flavonoid derivatives ( 1 ̶ 14 ) adhere to Lipinski, Veber, Egan, Muegge, Ghose and guidelines.Flavonoids ( 14) showed 1 violation in all the druglikeness parameters, making them unsuitable options for drug candidates.Lipinski's guidelines require that a compound satisfy five criteria before it can be administered orally, including a molecular weight (MW) of less than 500 g/mol, an octanol-water partition coefficient (logP) below 5, no more than 5 hydrogen bond donors (HBDs), no more than 10 hydrogen bond acceptors (HBAs), and a topological polar surface area (TPSA) of no more than 140 Å 2 .Veber's guidelines further specify that the number of rotatable bonds (nrotb) must be less than 10, and TPSA must be less than or equal to 140 Å 2 (consistent with Lipinski's guidelines) for good drug bioavailability.The study utilized SwissADME to assess the compliance of the most biologically active compounds ( 1 ̶ 14 ) with these guidelines.The results showed that all compounds complied with the specified Lipinski and Veber guidelines (see Table 2 ).In addition, the drug-like score of 1 indicated that all compounds met the criteria for developing novel medications with good theoretical evidence, as they exhibited high agreement (less than 6) for drug compounds with MW less than 500 g/mol, MLOGP 4.15, and Log S (ESOL) within the established guidelines.
Table 2 presents the physicochemical and drug-likeness properties of 16 flavonoids.Notably, flavonoid ( 9 ) contains the highest number of hydrogen bond donors and acceptors, while flavonoid ( 14) has the most rotatable bonds, which can improve oral bioavailability by increasing compound flexibility.Molecular weight is crucial to a compound's ability to reach its target site effectively.All flavonoids meet the criterion of having a molecular weight no greater than 450-500 g/mol.Flavonoids ( 1, 6 , and 7 ) have acceptable  a HIA: human intestinal absorption (%); BBB: blood -brain barrier penetration; PPB: plasma protein binding;CYP3A4: Cytochrome P4503A4; CYP2C19: Cytochrome P4502C19; hERG: human ether-a-go-go-related genehERG inhibition potential (pIC 50 ), the potential risk for inhibitors ranges from 5.5 − 6.
b The values are determined using admetSAR and SwissADME.
To efficiently develop a new drug, it is important to assess its pharmacokinetic properties, which include absorption, distribution, metabolism, excretion, and toxicity (ADMET).The study utilized admetSAR to evaluate the ADMET characteristics of the 14 active flavonoids ( 1 ̶ 14 ).The results are presented in Table 3 , which includes the evaluation of seven critical ADMET properties, including human intestinal absorption (HIA), blood -brain barrier (BBB) penetration, water solubility, cytochrome P450 enzyme (CYP3A4 and CYP2C19) inhibition, and hERG inhibition.It is important to note that crossing the blood -brain barrier can be beneficial for drugs that target the central nervous system (CNS) but may be harmful to those with little effect on the CNS.BBB penetration is classified as high ( > 2), medium (2-0.1), or low ( < 0.1).Table 3 shows that all flavonoids, except flavonoid ( 14) , have high HIA values, with ( 1, 8 , and 9) having complete absorption rates.However, none of the compounds, except flavonoids ( 8) and ( 9) , can penetrate the BBB.Furthermore, all flavonoids have good water solubility.Regarding cytochrome P450 enzymes, all flavonoids except flavonoid ( 1) show positive/negative inhibition of CYP3A4, while all flavonoids show negative/positive inhibition of CYP2C19.In addition, all flavonoids, including the control, show negative inhibition of the hERG gene.

Analysis of bioactivity and pharmacological properties of flavonoids
The study also utilized Molinspiration Chemoinformatics to evaluate the interactions between flavonoids and various biological targets, including G protein-coupled receptors (GPCRs), ion channel modulators (ICMs), kinase inhibitors (KIs), nuclear receptor ligands (NRLs), protease inhibitors (PIs), and enzyme inhibitors (EIs).A bioactivity score of 0.00 indicates a molecule's potential for exhibiting good downstream biological activity.On the other hand, a score between − 0.50 and 0 suggests good biological activity, and a score less than − 0.50 indicates inactivity.Flavonoid ( 9) from Table 4 demonstrated the best activity for all parameters and is the most promising candidate.Except for flavonoids ( 11 ̶ 13 ), all flavonoids showed good bioactivity in at least three mechanisms, thus narrowing down the list of potential drugs to flavonoids ( 7 ), which also has an excellent score in Table 2 .Overall, the findings indicate that the pharmacological activity of drug compounds may involve multiple mechanisms.Their interactions with GPCR ligands, nuclear receptor ligands, and inhibitors of proteases and other enzymes may play a crucial role in their physiological effects.

Molecular docking and postdocking analysis
The process of molecular docking is critical in determining the binding affinity between a ligand and a protein or receptor [93][94][95] .This method provides valuable insights into the behavior of compounds and their interactions with the active site residues of a protein, as well as downstream cellular processes.The docking results for 14 flavonoids with six proteins are presented in Table 5 , with HER2 (PDB: 7JXH) showing higher binding affinity with all flavonoids than the other proteins.A recent study carried out by Sajjan et al. showed HER2 (PDB: 7JXH) binding affinity to be − 8.5 kcal mol − 1 ), which is less than our obtained value [90][91][92][93][94][95] .
Flavonoid ( 1) displayed a high binding affinity with EGFR, HER2, FPPS, HGDS, and Keap1 on the Nrf2 protein.At the same time, flavonoid (6) demonstrated higher binding affinity than the control with EGFR, HER2, FPPS, and the highest affinity of − 9.4 kcal mol − 1 with (PDB:1 × 2 J) Keap1 protein on Nrf2.The results of (PDB:1 × 2 J) Keap1 protein on Nrf2 were better than the experimental study carried out by Cai et al., where they obtained a value of (− 9.1 kcal mol − 1 ) [93][94][95] .Flavonoids ( 8) and ( 12) had a lower affinity with four out of six proteins, possibly because they have halogens attached to the phenyl ring.The molecular docking study revealed The values obtained for all flavonoids for Nrf2 (PDB:1 × 2 J) were similar to those obtained for the complex compound with DHPA by Cai and his team [90][91][92] and the results are given in Table 5 .Researchers can optimize compounds and develop safe and effective therapeutic agents by combining computational and experimental methods.
In addition, analyzing the chemical bonds and interactions between the ligand and protein can help explain how the ligand affects the active parts of disease-causing pathogens.The study provides information on bond lengths for each type of bond and residue number.Figs. 3 and 4 show 2D diagrams for flavonoid ( 7) interaction with HPGDS (PDB: 1V40), DCK (PDB: 1P60) and Nrf2 (PDB: 1 × 2 J), including (a) the ligand in the protein pocket, (b) active site residues, (c) hydrogen bond density, and (d) ligand -protein interaction.The remaining docking pictures can be found in the supplementary information (see Figs. S15-S16 in the SI).The postdocking analysis of all 14 flavonoids and their respective proteins was conducted using various software tools, including UCSF Chimera, PyMOL version 2.5, and BIOVIA Discovery Studio.In Fig. 3 , the interacting active site residues are highlighted in orange, while heteroatoms are color-coded, with nitrogen in blue, hydrogen in white, and oxygen in red.  of the active site residues and their interactions with the ligand.The active site amino acids that interact with the ligand include Asp, Ser, Thr, Leu, Glu, Lys, Ile, Val, and Ala (for details, refer to Table 6 and S16 in the Supplementary Information).In this binding, a total of one hydrogen bond was identified between the ligand and the (PDB:1 × 2 J) active site, with one of them being particularly strong, with a distance of 2.207 Å. Fig. 3 c displays various intramolecular bonds in different colors, indicating interactions with the ligand.Two conventional hydrogen bonds occur between Ser and Asp with the oxygen of the phenyl ring, and eight Pi-alkyl bonds form between Ala, Ile, Leu, and Lys on the phenyl ring.Additionally, Thr is involved in two van der Waals forces with two oxygen atoms of the phenyl ring.Fig. 3 d highlights the hydrogen bond density within the ligand -protein pocket, with magenta representing donor hydrogen bonds and parrot green representing acceptor hydrogen bonds, as shown in Fig. 4 .In Fig. 4 , depicts the binding pocket of DCK (PDB: 1P60) containing flavonoid ( 7 ), while Fig. 4 b showcases the active site residues interacting with ( 7 ).The active site amino acids interacting with the ligand comprise Trp, Leu, Gly, Tyr, Cys, Met, Arg, and Asp.In this binding, seven hydrogen bonds were detected between the ligand and the (PDB:1P60) active site.Fig. 4 c uses various colors to represent different intramolecular bonds, including van der Waals, hydrogen carbon bonds, Pi-alkyl, and conventional hydrogen bonds, which interact with the ligand.Fig. 4 d displays the hydrogen bond density within the ligand -protein pocket.

Molecular dynamics (MD) simulation
To understand the stability and interactions of the protein -ligand complex of flavonoid ( 7) with Nrf2 (PDB:1 × 2 J) and DCK (PDB:1P60), molecular dynamics simulations were performed for 100 ns, as shown in Fig. 5 .The MD simulation results were analyzed using several metrics, such as RMSD, RMSF, Rg values, potential energies, temperature, and hydrogen bonding.The RMSD values of the protein -ligand complex, ligand, and protein backbone were calculated during the molecular dynamics simulation, as shown in Fig. 5 .The results indicated that the RMSD values for (PDB: 1P60), including the protein -ligand complex, ligand, and water ions, ranged between 0.05 and 0.4 nm ( Fig. 5 a-c).Meanwhile, the RMSD value for (PDB:1 × 2 J) was within 0.04 to 0.21 nm ( Fig. 5 a-c), suggesting that the ligand -protein complex underwent significant structural changes or movement, possibly due to the ligand binding event or conformational change induced by ligand -protein interactions.More detailed information can be found in the supplementary  document (see Figs. S17-S18 in the SI).Fig. 5 c displays the RMSD graph of the combined proteins, with the red curve representing the fluctuations of (PDB:1P60) and the black curve depicting the fluctuations of (PDB:1 × 2 J).Undoubtedly, (PDB:1P60) has a higher fluctuation and a steep increase at ∼45,000 ps, while (PDB:1 × 2 J) has no steep increase in fluctuation.This means that both proteins have different coordinates and different structural conformations, indicating differences in stability [92][93][94] .
The RMSF measures the average deviation of individual atoms in a protein from their average positions.In this simulation, separate RMSF values were calculated for the protein and the ligand.Fig. 6 displays the RMSF values of the protein -ligand complex of Nrf2 (PDB:X2J) and DCK (PDB:1P60) with flavonoids ( 7 ).The RMSF value for the complex (PDB: 1P60) and flavonoid ( 7) is approximately 0.04 nm with a slight fluctuation in the 700-800 range amino acid atom having a value of approximately 0.14 nm.These findings suggest that the protein and the ligand are relatively stable during the simulation, with only minor fluctuations in atom position.However, fluctuations in certain ligand atoms could be significant and potentially linked to a specific function or interaction of the protein.For the complex (PDB: 1 × 2 J) and flavonoid ( 7 ), the RMSF value ranges from 0.004 to 0.055 nm, indicating massive fluctuations in certain regions of the protein or the ligand, while others remain relatively rigid.High RMSF values could suggest conformational changes, flexibility, or disorder in the protein or the ligand.These alterations may be due to various factors, such as ligand binding, solvent exposure, or changes in temperature or pH [94][95][96] .The radius of gyration (Rg) in molecular dynamics is used to monitor changes in the shape and compactness of a molecule over time or to compare the shapes of different molecules.The Rg values for the ligand -protein complex between flavonoid ( 7) and (PDB:1P60) suggest an extended conformation, while the Rg values for the ligand alone suggest a compact conformation (see Fig. 7 a).The Rg values for the ligand -protein complex between flavonoid ( 7) and (PDB: 1 × 2 J) indicate a large and extended complex (see Fig. 7 ).According to the MD simulation results presented in ( Fig. 8 ) , the number of hydrogen bonds formed between compound ( 7) and the active sites of (PDB: 1P60) and (PDB: 1 × 2 J) ranged from 0 to 2 during the 100 ns simulation period.This variability suggests that the interaction between the ligand and protein is dynamic and that the ligand can adopt different conformations and binding modes within the protein-binding pocket.The formation of stable hydrogen bonds indicates that flavonoid ( 7) has the potential to interact with both proteins' active sites.
The temperature curves plotted in ( Fig. 9 ) demonstrated stability, ranging from 298.5 to 301.5 Kelvin throughout the 3000 PS simulation.This suggests that the system maintained consistent thermal energy, despite the presence of the protein -ligand complex.The potential energy of the system showed fluctuations ranging from − 5.1e + 05 to − 5.04e + 05 kJ mol ̶ 1 , indicating changes in electrostatic interactions between atoms within the system.
The simulation results suggest that compound ( 7) has the potential to act as an inhibitor of both protein DCK and Nrf2, as it showed a stronger interaction with (PDB:1P60) and (PDB:1 × 2 J).This finding is promising, given that mRNA expression of DCK has been associated with poor prognosis in liver cancer patients.Therefore, flavonoid ( 7) could potentially be a therapeutic agent for liver diseases.However, further experimental studies are necessary to verify its efficacy and safety as a potential liver therapeutic.
To validate our findings and determine whether flavonoid ( 7) is optimal, we extended our investigation to include two additional flavonoids ( 10 and 14 ) and the standard drug, Taxifolin, which exhibited the most favorable docking results.Subsequently, we conducted molecular dynamics simulations, and the outcomes are presented in Fig. 10 , comprising three panels: a) ligand comparison, b) protein comparisons, and c) protein-complex comparisons.Fig. 10 shows the root mean square deviation (RMSD) analysis for the three flavonoids and the standard drug.In Fig. 10 c the red line corresponds to flavonoid (7) , while the blue, green, and black lines represent flavonoids ( 10, 14) , and the standard drug, respectively.Flavonoid (7) is the only one exhibiting an RMSD value close to 5, whereas the others collectively yield a value of 6.4.This notably higher combined value indicates that these ligands possess dissimilar shape orientations and exhibit less structural similarity.Notably, flavonoid ( 7) displays a lower RMSD value than the other compounds, as evident in Figs. 10 a and 10 b.This underscores its unique structural characteristics and suggests its potential as the most promising candidate.For more detailed information, including individual graphs for each flavonoid with respect to the reference molecule (PDB:1 × 2 J), please refer to the supplementary information (Figs.S17-S21).
Fig. 11 illustrates the merged radius of gyration for three flavonoids and the standard drug.In Fig. 11 (a) we observe a radius of gyration range of 0.34-0.35nm for flavonoid 7 , 0.37-0.38nm for the standard drug, and 0.40-0.42nm for flavonoids ( 10 and 14) .It is important to note that higher values in the radius of gyration indicate lower structural similarity and fewer interactions among molecules.Remarkably, flavonoid ( 7) exhibits the highest structural similarity and interactions among its molecules due to its lower radius of gyration.Detailed radius of gyration values and graphs for each flavonoid can be found in the supplementary information (see Figs. S22-S24).
Fig. 12 depicts the root mean square fluctuation (RMSF) values for four flavonoids, including the standard drug.RMSF serves as a metric for assessing the flexibility of a structure by evaluating the squared deviations from a reference structure.It is important to note that RMSF values are context-dependent and influenced by factors such as temperature, pressure, solvent conditions, and the choice of reference structure, which can vary between an average, crystallographic, or initial structure.These contextual and reference-related factors can lead to variations in RMSF values for the same structure, underscoring the need to consider these factors carefully.In general, an RMSF value within the range of 0.05 to 0.57 nm signifies varying degrees of flexibility, spanning from low to high, compared to other structures exposed to similar conditions.In our specific analysis, all RMSF values fall within this designated range.Notably, flavonoids ( 7) and ( 14) exhibit a peak at 0.57 nm, whereas flavonoid 10 demonstrates a peak at 0.54 nm, and the standard drug exhibits the lowest peak at 0.45 nm.They acknowledge that the interpretation of the RMSF values for flavonoids ( 7) and ( 14) hinges on their comparison with other structures within the same context, and the reference framework is crucial.For instance, the standard drug's RMSF value of 0.05 nm sets a benchmark for low flexibility.Consequently, structures with RMSF values ranging from 0.5 to 0.6 nm are highly flexible.For a more comprehensive understanding of the RMSF values for each protein complex, please refer to the supplementary information provided in Fig. S25 .
PCA was conducted on the MD trajectories of flavonoid ( 7 ) with protein (PDB:1P60) and the (PDB:1 × 2 J) protein -ligand complex at 300 K using the Bio3D package, as shown in Fig. 13 .The results in Fig. 13 show that PC1 captured the highest variability for (PDB:1P60) (21.29%) and (PDB:1 × 2 J) (12.04%), indicating significant internal motions within the MD trajectories.PC2 had lower  variance percentages (6.91%) for 1P60 and (9.07%) for the (PDB:1 × 2 J) protein complex, while PC3 also explained 21.29% and 6.91% of the variance, respectively.Together, these three components explained a cumulative variance of 49.49% for the 1P60 complex and 33.15% for the (PDB:1 × 2 J) complex.These findings offer valuable insights into the collective motions and conformational changes occurring in protein -ligand complexes during MD simulations.Additionally, the cosine content values of the eigenvectors calculated from the MD trajectory for the (PDB:1 × 2 J) and (PDB:1P60) complexes were determined to be 0.59 and 0.66, respectively.These values indicate that the simulations have reached convergence, demonstrating the stability and reliable sampling of the protein -ligand complexes throughout the MD simulations.Furthermore, PCA was conducted on flavonoids ( 10 and 14) , and the standard drug (taxifolin) in conjunction with the protein (PDB: 1 × 2 J) is shown in Fig. 14 .This additional analysis was performed to validate our study's results and ascertain whether flavonoid ( 7) is a suitable candidate for liver drug development.It is worth noting that relying solely on docking results may not conclusively establish flavonoid ( 7) as the optimal choice.Consequently, flavonoids ( 10 and 14) and the standard drug were selected to examine their molecular dynamics trajectories.The cosine content for flavonoids ( 10 and 14) was calculated to be 0.01, while for the standard drug (taxifolin), it was 0.70.Cosine content values approaching 0 suggest that the molecular simulation explores diverse regions within the conformational space, indicating complex or irregular motion.However, it is crucial to recognize that these interpretations are not absolute and should be considered alongside other factors, including the simulation duration, system size, and choice of coordinates.Moreover, the cumulative variance of the three principal component analysis (PCA) components for flavonoids ( 10) and ( 14) in conjunction with protein (PDB:1 × 2 J) amounts to 39.02%.This finding implies that flavonoid ( 7) exhibits a lower cumulative variance than the others, indicating a more simplified or coarse-grained trajectory representation, as shown in Fig. 14 .

Conclusion
This study investigated analogs of flavones ( 1 -14 ) to determine their efficacy as liver drugs.These derivatives ( 1 -14 ) have been investigated computationally for their potential use against three cancer-related proteins and one parasite target.This work includes DFT calculations, molecular docking calculations, binding energy calculations, thermodynamic properties, HOMO and LUMO investigation, drug-likeness analyses, ADMET properties, MD simulation, and their biological activity spectra.Among all fifteen compounds, flavonoid ( 7 ) showed good agreement binding affinity as a liver therapeutic, as it has shown excellent efficiency through computational analysis.These results indicated that 5,7-dihydroxyflavone might be a promising drug molecule, as it showed good binding affinity compared to the standard drug 7,8-dihydroxyflavone with EGFR, HER2, FPPS, HGDS, and Keap1 on the Nrf2 protein.Moreover, flavonoids ( 7 ) obeyed the Lipinksi, Ghose, Veber, Egan, and Muegge rules.Flavonoid ( 7) showed an excellent intestinal absorption rate.However, validating all the results in the wet laboratory would be the best way to characterize the properties of these compounds as therapeutic liver drugs.

Declaration of Competing Interest
The author declares that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2
displays the HOMO-LUMO orbitals of flavonoids ( 7 ), while the remaining figures can be found in the Supporting Information (see Figs. S1 to S5).NBO charge structures for compounds ( 1 ̶ 14 ) are shown in Figs.S6 to S11 in the SI.

Fig. 5 .
Fig. 5.The RMSD evolution (a) for the docked complex, green line (flavonoid 7 , red line with PDB: 1P60, black line); (b) for the docked complex, green line (flavonoid 7 , red line with PDB: 1 × 2 J, black line); and (c) for the merged docked complex between two proteins (PDB: 1P60, red line and PDB: 1 × 2 J, black line) during the 100 ns MD simulation.

Fig. 8 .
Fig. 8.The number of hydrogen bonds versus time (ps) plots for hydrogen bond stabilization hydrogen bonding for (a) the protein complex of PDB: 1P60 and flavonoid 7 and (b) the protein complex of PDB: 1 × 2 J and flavonoid 7 during a 100 ns MD simulation.

Fig. 9 .
Fig. 9.The temperature and potential energy curves over the course of the 100 ns MD simulation for (a) a potential energy graph for the protein -ligand complex of PDB: 1 × 2 J and flavonoid 7 and (b) a temperature graph of the protein -ligand complex of PDB: 1 × 2 J and flavonoid 7 .(c) Potential energy graph for the protein -ligand complex of PDB: 1P60 and flavonoid 7 .(d) Temperature graph of the protein -ligand complex of PDB: 1P60 and flavonoid 7.

Fig. 7 a
shows the radius of gyration for protein complex (PDB:1P60) with flavonoid ( 7) , and Fig. 7 b shows the radius of gyration values of flavonoid 7 with (PDB:1 × 2 J) protein.Fig. 7 c shows Rg values of two proteins merged (PDB:1 × 2 J and PDB:1P60) with ligand, water ion and ligand -protein complex; in both lines, the straight line at 0.3 nm indicates water ions.

Fig. 12 .
Fig. 12. MD simulation evolution of RMSD for merged graphs of docked protein -ligand complexes between PDB proteins: 1 × 2 J and modeled ligands 2 (blue line), 10 (red line), and 14 (black line) and the reference drug (taxifolin) (green line) during 100 ns of MD simulation.

Fig. 13 .
Fig. 13.PCA conducted on the MD trajectory of the (a) 1P60 protein and (b) PDB:1 × 2 J protein at 300 K.In the plot, the blue dots correspond to energetically unstable conformational states, while the red dots indicate stable conformational states.

Table 2
In silico prediction of physicochemical parameters for the flavone and its derivatives ( 1 -14 ). a HBA, number of hydrogen bond acceptors; nroth, number of rotatable bonds; TPSA, topological polar surface area (Å 2 ); standard reference drugs such as tropoflavin and taxifolin.

Table 3
In silico prediction of selected ADMET parameters for all flavones ( 1 -14 ). a

Table 5
Docking results for flavone derivatives ( 1 ̶ 14 ) against six protein targets.a